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I.  Introduction 


The  Modal  Superposition  Method  is  a  useful  algorithm  for  the 
solution  of  linear  problems  in  structural  dynamics.  The  basic  idea  of 
the  method  is  to  transform  the  coupled  equations  of  motion  into  an 
uncoupled  form.  This  change  of  basis  is  equivalent  to  diagonalization 
of  the  underlying  matrix  differential  equation,  and  has  the  effect  of 
trading  one  complicated  n-dimensional  system  of  ODE  for  n  simple 
(scalar)  ODE.  Typically,  this  change  of  basis  to  modal  coordinates  is 
very  efficient  computationally,  and  linear  problems  are  generally 
solved  in  this  manner.  Often,  not  all  the  modes  participate  in  the 
solution,  and  so  the  problem  may  be  projected  onto  a  smaller  subspace, 
with  a  corresponding  decrease  in  computational  cost. 

The  principle  of  superposition  does  not  apply  to  nonlinear 
problems,  so  Modal  Superposition  is  not  strictly  valid  in  a  nonlinear 
setting.  The  method  can  be  extended  to  this  more  difficult  case  by 
considering  the  linearization  of  the  nonlinear  problem  in  a  local 
region.  This  linearized  problem  can  be  used  to  generate  modal 
coordinates,  and  these  modes  can  be  used  as  generalized  coordinates  in 
the  solution  of  the  nonlinear  problem.  The  presence  of  nonlinearities 
has  the  effect  of  introducing  coupling  among  the  modal  equations  of 
motion  but  the  reduction  in  size  of  the  problem  obtained  by  projection 
to  a  smaller  modal  basis  typically  justifies  the  increased  cost  due  to 
this  coupling.  In  addition,  the  knowledge  that  off-diagonal  terms 
will  exist  in  the  modal  problem  means  that  the  exact  modes  (the  ones 
that  diagonalize  the  linearized  problem)  need  not  be  found:  an 
approximation  will  do,  since  the  effect  of  the  error  is  to  add 
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off-diagonal  terms  to  the  (already  ooupled)  modal  problem. 

The  scope  of  this  report  is  to  propose  a  general  algorithm  for  the 
solution  of  a  class  of  nonlinear  problems  in  structural  dynamics  using 
a  modal  representation.  The  development  of  some  of  the  underlying 
linear  theory  will  be  presented!  followed  by  generalizations  to 
nonlinear  problems.  These  topics  will  be  used  to  develop  an  algorithm 
for  the  nonlinear  case,  and  several  one-dimensional  examples  will  be 
presented.  This  research  is  a  starting  point  in  an  attempt  to  produce 
an  adaptive  algorithm  for  the  solution  of  nonlinear  problems,  and  so  a 
number  of  topics  for  further  development  will  be  presented. 

It  is  assumed  that  the  reader  has  a  basic  knowledge  of  structural 
dynamics,  finite  element  methods,  and  matrix  spectral  theory. 
Sufficient  detail  on  these  topics  can  be  found  in  references  [1],  [2], 
and  [3]. 

Some  consideration  of  semantics  should  be  noted:  in  the  following 
development,  a  "mode”  is  taken  to  be  a  generalized  coordinate  vector 
that  reflects  some  aspect  of  the  vibrational  behavior  of  the 
structure,  but  not  neoessarlly  an  exact  eigenvector  of  the 
mathematical  model.  (This  is  an  ambiguous  definition,  but  the  word 
"mode"  is  considerably  more  imprecise  than  "eigenvector" ,  with  which 
it  is  often  taken  to  be  synonymous . }  Different  modes  are  taken  to 
satisfy  some  type  of  orthogonality  requirement  (usually  with  respect 
to  the  finite  element  mass  matrix),  and  to  have  unit  length  in  some 
appropriate  norm,  but  they  do  not  have  to  produce  a  diagonal  form  of 
the  equations  of  motions,  as  in  linear  theory.  This  convention  is 
merely  for  the  convenience  of  lumping  classes  of  generalized 
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coordinates  into  the  term  "mode",  instead  of  referring  to  "exact 
eigenvectors,  Ritz  vectors,  Lanczos  approximations  to  eigenvectors", 
eto. 


Finally,  all  important  symbols  are  tabulated  in  Appendix  1. 


II.  Linear  Modal  Analysis 


The  equations  of  motion  for  a  linear-elastic  structure  are  a  set  of 
partial  differential  equations  involving  both  spatial  and  temporal 
derivatives.  These  PDE  can  be  cast  in  a  semi-analytic  form  by 
performing  a  finite  element  discretization  in  the  space  domain.  This 
process  turns  the  spatial  differential  terms  into  an  algebraic  set  of 
equations.  The  time  derivatives  are  left  in  continuous  form,  and  the 
resulting  system  of  ordinary  differential  equations  can  be  written  for 
an  undamped  system  as: 

Kx  ♦  Mx  s  f  Eqn  1 

where  K  and  M  are  the  finite  element  stiffness  and  mass  matrices,  x  is 
the  veotor  of  nodal  displacements,  arid  f  is  a  load  vector  which  is  a 
function  of  time.  The  modal  coordinates  of  this  system  can  be  found 
by  solving  the  generalized  eigenvalue  problem: 


(K  -  w2M)u  *  0 


Eqn  2 
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The  pair  (K,M)  is  termed  a  matrix  pencil  (this  term  arises  from  the 
study  of  geometric  optics).  In  general,  K  and  H  are  positive-definite 
matrices,  and  so  the  eigenvalues  v?  are  real  and  positive.  Under 
these  conditions,  a  complete  set  of  eigenvectors  exists  that  are 
orthonormal  with  respect  to  the  mass  matrix: 

U  s  [u1,  U2,  •••  un]  Eqn  3 
UtMU  s  I  UtKU  s  D  r  diag(w,2  ) 


The  change  of  basis  to  modal  coordinates  z  is  given  by 


Dz  i  x  Eqn  4 


These  equations  can  easily  be  integrated  to  give  the  modal 
coordinates  z,  which  are  then  transformed  using  Eqn  4  to  give  the 
desired  solution.  Any  numerical  integration  scheme  appropriate  for  a 
single-degree-o f-freedom  system,  such  as  the  Duhamel  Integral  (see 
[1],  Pgs  100-108)  can  be  used  to  solve  the  modal  equations. 
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In  practice,  the  solution  of  linear  problems  generally  includes  a 
damping  term  Ci  to  account  for  energy  loss  in  the  structure.  In  order 
to  diagonalize  K,  C,  and  M  simultaneously,  the  damping  matrix  C  is 
generally  taken  to  be  some  combination  of  powers  of  K  and  M.  (In  many 
nonlinear  problems,  which  are  the  ultimate  scope  of  this  research,  the 
stiffness  K(x,  t)  typically  includes  nonconservative  terms,  and  so  the 
construction  of  a  damping  matrix  is  not  considered  further. )  The 
presence  of  damping  in  a  linear  system  causes  the  transient 
free-vibration  response  of  the  system  (corresponding  to  the  homogenous 
solution  of  Eqn  1)  to  decay,  leaving  the  steady-state  forced  response 
induced  by  the  load  f  ( the  nonhomogenou?  solution  of  Eqn  1 ) .  In  the 
following  discussion,  the  "response"  of  the  structure  is  taken  to  be 
the  steady-state  response. 

The  loading  term  ujt  represents  a  generalized  load  acting  on  the 
ith  modal  coordinate.  If  U|  and  f  are  nearly  orthogonal,  then  this 
load  is  small  relative  to  other  generalized  modal  loads.  In  this 
case,  the  ith  mode  does  not  appreciably  participate  in  the  solution, 
and  can  be  neglected  in  the  analysis.  Typically,  in  the  modelling  of 
a  complex  structure,  most  modes  do  not  participate,  and  thus  a  reduced 
problem  can  be  solved  by  projecting  the  solution  onto  a  p-dimensional 
subspace  spanned  by  the  p  "most  important"  modes: 

x  s  Upz  Eqn  6 


In  this  case,  the  columns  of  Up  (an  nxp  matrix)  are  the  p  dominant 
orthonormal  modes,  and  z  (a  pxl  vector)  is  the  reduced  set  of  modal 
coordinates.  It  is  important  to  observe  that  U^Up  >  Ip,  the  identity 
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matrix  for  p-dimensional  space,  but  that  UpOp  f  In,  unless  p  =  n.  In 
spite  of  this  last  fact,  Dp  is  still  often  termed  an  orthogonal 
matrix,  and  this  terminology  will  be  used  herein. 

The  load  vector  f  can  be  decomposed  into  orthogonal  components 
f  a  g  +  h  where  g  lies  in  the  column  space  (or  range)  of  Dp,  and 
gTh  s  0.  The  vectors  g  and  h  are  defined  by  projecting  f  onto  the 
range  of  Dp: 

g  =  OpOjf  h  =  (I  -  DpOp ) f  Eqn  7 

The  Pythagorean  theorem  implies  that  gTg  +  hTh  =  fTf  (i.e. 
||g||2  ♦  ||h||2  s  ||f||2)»  so  that  it  is  to  be  hoped  that  gj&  »  hTh  in 
order  that  the  reduced  set  of  coordinates  gives  a  good  approximation 
to  the  load.  Finally,  in  preparation  for  subsequent  development,  it 
should  be  noted  that  the  projection  defined  by  Eqns  6  and  7  is  valid 
for  any  matrix  Dp  with  orthonormal  columns,  but  that  the  particular  Dp 
formed  using  the  eigenvectors  of  the  pencil  (K,  M)  diagonalizes  the 
linear  problem. 

Once  the  important  modes  of  the  linear  problem  have  been  found,  the 
scalar  ODE  of  Eqn  5  oan  be  integrated  to  give  the  solution.  A  problem 
oocurs  when  the  time  scales  defined  by  the  modal  frequencies  v\  span 
several  orders  of  magnitude: 

max(ws )  »  min(w> ) 
i  1  i  1 

In  this  oase,  any  step  size  appropriate  for  numerical  integration  of  a 


"fast"  problem  (max(w:))  will  be  too  small  for  efficient  integration 

■  ® 


problem  is  inaccurate  for  a  fast  one.  If  the  load  function  f  defines 
smother  time  scale  (i.e.  the  duration  or  the  rapid  variation  of  an 
earthquake  record),  then  the  situation  is  even  more  complicated.  In 
any  case,  a  wide  range  of  time  scales  exists  in  the  problem,  and  any 
computational  scheme  must  accomodate  this  difficulty.  Differential 
equations  characterized  by  this  range  of  scales  are  referred  to  as 
stiff:  in  this  case,  the  stiffness  is  inherent  in  the  spread  of  the 
eigenvalues  w*  of  the  pencil  (K,M).  In  a  linear  modal  analysis,  the 
stiff  nature  of  the  equations  requires  selection  of  an  appropriate 
step  size  for  integration  of  the  modal  equations:  since  these 
equations  are  uncoupled,  these  step  sizes  may  be  chosen  independently. 
An  algorithm  that  does  not  completely  decouple  the  modes  may  require 
short  steps  for  a  slow  mode  if  it  remains  coupled  to  a  faster  one. 
Thus,  there  is  strong  incentive  to  diagonalize  (or  nearly-diagonalize ) 
any  linear  problem,  except  in  very  special  circumstanoes  (such  as  an 
extremely  short-duration  loading  and  response,  where  the  cost  of 
finding  the  modes  cannot  be  amortized  over  sufficient  time  steps  to 
remain  competitive). 

Two  algorithms  for  generation  and  solution  of  reduced  equations  for 
linear  problems  have  appeared  recently  in  the  literature.  Although 
these  algorithms  are  motivated  from  different  basic  principles,  in 
fact  they  are  closely  related.  Each  represents  an  efficient 
alternative  to  solution  of  the  problem  using  exact  eigenvectors. 
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The  first  approach  is  due  to  Wilson,  et.  al.  [4] [5],  who  suggests 
using  so-called  "Ritz  vectors"  to  achieve  a  reduced  set  of 
coordinates.  The  identification  of  these  Ritz  vectors  proceeds  from  a 
physical  argument  involving  the  solution  of  a  static  problem,  followed 
by  corrections  accounting  for  dynamic  terms  missed  in  previous  steps. 
(The  reader  is  referred  to  [4]  for  the  details.  The  following  summary 
involves  some  notational  changes  from  this  original  reference. )  In 
summary,  the  Ritz  vectors  q;  are  obtained  using  the  recurrence: 


1) 

Kx-) 

II 

II 

2) 

=  x-(/(xJmx-|  ) 

3) 

For 

i  a  2  to  p 

3a) 

KXj  =  Mxj.-j 

3b) 

*i  a  xi  “Sj<i(q[Mxi  )qi 

3c 

q,  =  yj/CyjMy.  ) 

(f  =  equivalent  static  load) 


i 


The  nxp  matrix  Qp  =  [q1iq2I  ...  Iqp]  is  (in  the  absence  of 
round-off  error)  an  orthogonad  matrix.  In  practice,  step  3b  of  the 
algorithm  is  badly  affected  by  round-off  error  in  a  low-precision 
computationad  environment,  and  needs  to  be  replaced  with  a  more  robust 
orthogonalization  scheme  (such  as  the  one  found  in  [33*  pg  105-110), 
in  order  to  achieve  orthogonality  of  the  columns  of  Qp  to  working 
precision.  If  enough  vectors  are  included  in  Qp,  these  reduced 
coordinates  give  excellent  results:  in  the  example  problems  studied 
in  [4],  this  algorithm  was  clearly  superior  to  the  use  of  exact 
eigenvectors.  The  authors  suggest  that  further  savings  could  be 


realized  by  diagonalizing  the  reduced  stiffness  QpRQp  before  numerical 
integration  is  performed.  Since  p  «  n,  this  last  step  is  very 
inexpensive. 

Examination  of  this  algorithm  reveals  that  the  Ritz  vectors  are 
chosen  to  be  an  orthonormal  basis  for  the  space  spanned  by  the  set 
UrV,  (K'1M)K‘1f,  (rVrV,  ...  ,  (^Mf'VfJ.  (This  subspace  is 
termed  the  Krylov  subspace  of  order  p,  generated  by  the  matrix  K^M  and 
the  vector  K*V:  reference  [3]  contains  several  excellent  chapters  on 
the  utility  of  these  Krylov  subspaces.)  The  algorithm  of  references 
[4]  and  [5]  therefore  generate  these  Ritz  vectors  by  direct 
calculation  of  an  orthonormal  basis  for  this  subspace,  with  full 
orthogonalization  (step  3b)  at  each  step. 

The  second  approach  for  generation  of  a  reduced  problem  is  due  to 
Nour-Omid  and  Clough  [6].  This  method  uses  the  Lanczos  Algorithm 
(refer  to  [3]  or  to  [7]  for  a  detailed  derivation)  to  generate  Lanczos 
veotors:  tnese  vectors  define  the  reduced  coordinate  set  for  the 
problem.  The  Lanczos  Algorithm  involves  the  use  of  a  three-term 
recurrence  to  generate  an  orthonormal  basis  for  the  Krylov  subspace 
mentioned  in  the  last  paragraph.  This  recurrence  (even  after 
modifications  are  made  to  reduce  the  effect  of  round-off  error)  is 
considerably  less  expensive  than  the  orthogonalization  scheme  used  in 
C4] [5],  and  results  in  a  tridiagonal  form  for  the  equations  of  motion. 
The  resulting  reduced  problem  can  easily  be  diagonalized,  or  a 
tridiagonal  formulation  for  the  numerical  integration  of  the  equations 
of  motion  can  be  used.  In  either  oase,  the  computational  cost  is  much 
less  than  any  method  for  solving  the  full  (unreduced)  problem. 
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III.  Nonlinear  Modal  Analysis 


A  common  class  of  nonlinear  problems  in  structured  dynamics  occurs 
when  the  stiffness  and/or  mass  matrices  are  functions  of  the 
displacement  x  (or  its  temporal  derivatives).  In  this  case,  the 
equilibrium  equations  given  by  Eqn  1  become: 

K(x,x,jc,t)x  +  M(x,x,x,t)x  =  f(x,x,x,t)  Eqn  8 


The  most  typical  example  of  this  class  of  problems  is  considerably 
simpler: 

K(x,t)x  Mx  =  f(x,t)  Eqn  9 

In  this  case,  the  mass  matrix  is  constant,  and  the  stiffness  and  force 
depend  only  on  x  and  t.  The  following  development  assumes  that  the 
nonlinearities  can  be  cast  in  the  form  of  Eqn  9:  the  more  general 
case  given  by  Eqn  8  can  be  analyzed  in  a  similar  manner. 

For  a  given  configuration  xq  =  x(t0),  the  pencil  (K(xQ,t0)  M)  can 
be  used  to  generate  a  set  of  modes.  These  modes  could  be  the  exact 
eigenvectors  that  diagonalize  this  pencil  or,  as  in  the  last  section, 
could  be  Lanczos  vectors  that  tridiagonallze  mtV  Since  the  problem 
is  nonlinear,  it  is  not  expected  that  this  diagonal  (or  tridiagonal) 
structure  will  be  preserved  for  x  i  Xq.  (To  be  more  precise,  the 
nonlinear  problem  can  be  linearized  at  xq,  and  a  set  of  reduced  modal 
coordinates  extracted. ) 
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The  primary  questions  in  such  a  local  reduction  are  the  number  of 
reduced  coordinates  required  to  gain  an  acceptable  approximation  in  a 
region  near  xq,  and  the  size  of  the  region  over  which  the 
approximation  is  valid.  Assuming  that  a  representative  load  vector 
(or  some  generalization  of  a  participation  factor)  can  be  found,  the 
methodology  of  references  [4],  15],  and  [6]  can  be  used  to  determine 
the  number  of  coordinates  required  at  xq.  The  magnitudes  of  the 
projected  load  and  the  error  involved  in  this  projection  can  be 
monitored  to  insure  that  enough  modes  have  been  included  to  maintain 
accuracy.  Some  measure  of  off-diagonal  strength  in  the  reduced 
equations,  such  as 


V(rji  r  j  j  ) 


(where  t*-,j  is  an  element  of  the  reduced  stiffness)  can  be  used  to 
determine  whether  the  modal  coupling  is  too  strong  (as  in  solving 
stiff  problems).  All  these  topics  are  currently  being  studied,  and 
the  results  will  be  presented  in  future  publications. 


For  some  problems  (see  reference  [8]  and  the  examples  presented  in 
the  following  sections),  good  results  have  been  obtained  by  taking  xq 
and  tg  to  be  zero  (i.e.,  using  the  initial  configuration  to  generate 
mode  shapes)  and  applying  these  same  reduced  coordinates  throughout 
the  entire  nonlinear  problem.  On  the  other  hand,  when  the 
nonlinearities  Involve  large  inelastic  effects,  it  is  expected  that 
the  modes  of  vibration  will  need  to  be  updated  at  some  stages  of  the 
solution  process. 


IV.  Proposed  Algorithm  for  Reduced  Coordinate  Analysis 


Using  the  direct  integration  of  the  equations  of  motion  as  a  model, 
and  including  some  scheme  (such  as  those  just  presented)  to  generate  a 
set  of  reduced  coordinates,  the  following  skeletal  algorithm  is 
presented  as  an  outline  of  a  general  solution  scheme  for  nonlinear 
problems  of  the  form  given  in  Eqn  9: 

1)  Initialize: 

la)  set  initial  conditions  on  x,  x,  t,  f 

lb)  calculate  initial  values  of  K,  it,  M 

lc)  evaluate  mode  shapes  Qp 

2)  For  t  a  to  to  t  s  t_jnax 

2a)  update  modes,  if  necessary 

2b)  form  reduced  problem  Rz  ♦  z  =  g 

2c)  integrate  reduced  problem  to  end  of  time  step 

2d)  if  equilibrium  satisfied  then  done  with  step 

else  apply  iteration  strategy  and  repeat  from  (2b)  or  (2a) 

This  algorithm  is  intentionally  nebulous,  sinoe  the  development  of 
this  scheme  is  still  proceeding.  As  mentioned  in  the  last  section, 
the  question  of  when  to  update  the  modes  is  the  subject  of  present 


research 
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V.  Computational  Costs  of  Reduced  Algorithm 


In  order  to  solve  the  equations  of  motion  (either  full  or  reduced), 
some  numerical  integration  scheme  such  as  Newmark's  Method  (see  [9]) 
is  generally  employed.  For  the  full  (unreduced)  set  of  coordinates, 
the  stiffness  matrix  must  be  factored  at  each  time  step:  this 
operation  is  of  the  order  n*b3,  where  n  is  the  number  of  equations  and 
b  is  an  average  bandwidth  of  the  stiffness  matrix.  For  a  reduced 
coordinate  system,  the  stiffness  matrix  must  be  factored  each  time  the 
coordinates  are  updated,  but  it  is  assumed  that  this  occurs 
infrequently,  compared  to  the  number  of  time  steps.  The  reduced 
stiffness  has  to  be  factored  at  each  time  step:  this  requires  on  the 
order  of  p3  operations,  where  p  is  the  number  of  reduced  coordinates. 
Formation  of  the  reduced  stiffness  using  the  least  efficient  means 
(i.e.,  forming  the  unreduced  stiffness  and  pre-  and  post-multiplying 
by  the  orthogonal  modal  matrix)  requires  on  the  order  of  n*b*p 
operations.  Therefore,  the  most  inefficient  reduced  formulation  will 
require  on  the  order  of  n*b#p  ♦  p3  operations.  If  p  is  much  less  than 
b  (which  is  typical  of  large  problems,  but  definitely  not  the  case  for 
the  one-dimensional  examples  presented  below),  then  the  reduced 
algorithm  will  be  competitive  with  direct  integration.  If  measures 
are  taken  to  improve  the  computational  efficiency  of  the  reduced 
algorithm  (see  [10]  for  one  example),  these  asymptotic  operation 
counts  for  reduced  analyses  should  improve  relative  to  the  direct 
formulation. 
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When  Iteration  within  each  time  step  is  considered,  the  comparison 
of  computational  cost  is  less  clear.  At  this  stage  of  the  research, 
it  is  premature  to  draw  conclusions  as  to  whether  the  reduction  of 
coordinates  enhances  or  retards  convergence  of  a  given  iteration 
scheme.  Development  of  a  robust  reduced  Iteration  scheme  for 
satisfaction  of  the  equations  of  equilibrium  at  the  end  of  each  time 
step  is  presently  one  of  the  main  topics  of  study  by  the  authors.  The 
small  order  of  the  reduced  equations  makes  oertain  iteration  schemes 
(such  as  those  in  the  Newton  family)  more  feasible. 


VI.  Example  Problems 


Several  one-dimensional  problems  are  analyzed  in  order  to  calibrate 
and  test  the  reduced  coordinate  models.  For  each  test  problem,  three 
types  of  analysis  are  performed: 

1)  Direct  integration  of  the  equations  of  motion 

2)  Integration  of  the  reduced  system, 

using  exact  eigenvectors  for  reduction 

3)  Integration  of  the  reduced  system, 

using  Lanczos  vectors  for  reduction 

In  the  second  and  third  cases,  the  reduced  coordinate  algorithm 
presented  in  Section  17  is  used.  In  all  cases,  the  equations  of 
motion  are  numerically  integrated  using  Newmark’s  Method  (see  [9]), 
with  integration  parameters: 


alpha  s  1/4 
delta  =  1/2 

The  first  example,  used  for  calibration  purposes,  is  that  of  free 
vibration  of  a  linear  string.  The  assumption  of  a  linear  displacement 
field  over  each  element  leads  to  an  initial  value  problem  in  the  form 
of  Eqn  1,  with  a  mass  matrix  equal  to  the  identity,  and  a  tridiagonal 
stiffness  matrix: 

K  s  tridiag(-1 ,  2,-1) 

The  tridiagonal  stiffness  is  equal  to  the  so-called  Jacobi  matrix: 
the  eigensystem  of  this  matrix  is  known  in  closed  form,  so  an  exact 
solution  is  available  for  comparison.  The  number  of  equations  chosen 
for  this  problem  is  fifty.  The  fundamental  period  of  this  system  is 
approximately  102,  and  the  numerical  schemes  were  integrated  from 
t  s  0  to  t  s  102,  using  time  steps  of  unit  length. 

The  first  ten  mode  shapes  are  shown  in  Figures  1  and  2.  Both  sets 
of  modes  were  obtained  from  the  Lanczos  algorithm:  the  exact 
eigenvectors  can  be  obtained  by  finding  the  eigensystem  of  the 
trldlagonal  matrix  that  gives  the  three-term  recurrence  used  by  the 
Lanczos  soheme.  Typically,  in  order  to  calculate  p  eigenvectors,  more 
than  p  steps  of  the  recurrence  must  be  taken,  so  finding  eigenvectors 
using  the  Lanczos  method  is  always  more  expensive  than  finding  an 
equal  number  of  Lanczos  vectors.  A  vector  with  each  component  set 
equal  to  unity  was  used  to  start  the  algorithm.  Note  that 
antisymmetric  modes  appear  in  both  cases.  Since  the  stiffness  and 
mass  are  persymmetric  (symmetric  about  both  diagonals),  if  a  symmetric 


starting  vector  is  used  to  generate  modes,  in  the  absence  of  round-off 
error,  only  symmetric  modes  would  be  produced.  Obviously,  one  effect 
of  round-off  error  is  to  include  antisymmetric  modes  that  do  not 
participate  in  the  solution. 

Because  the  problem  is  one-dimensional,  the  history  of  the 
displacement  can  be  visualized  in  three  dimensions:  space,  time  and 
displacement.  Figure  3  shows  the  solution  domain  in  the  space- time 
plane,  and  the  viewpoint  for  display  of  the  displacement  history.  One 
of  the  main  reasons  for  choosing  one-dimensional  problems  for  test 
purposes  (other  than  low  computational  cost)  is  that  the  entire 
history  of  the  problem  can  be  summarized  in  the  display  of 
displacement  as  a  function  of  space  and  time. 

Figures  4  thru  7  display  displacement  histories  for  the  various 
methods:  exact  solution,  direct  integration,  eigenvectors  for  reduced 
coordinates,  and  Lanczos  vectors  for  reduced  coordinates.  Figures  8 
and  9  show  the  displacement  errors  for  the  three  numerical  integration 
schemes,  including  results  using  six  and  ten  reduced  coordinates. 
These  errors  are  tabulated  in  Table  1 ,  along  with  estimates  of 
computational  costs.  Note  that  since  the  problem  is  tridiagonal,  the 
reduced  algorithm  is  placed  at  a  serious  disadvantage,  for  the  reasons 
given  in  the  section  on  computational  costs.  This  disadvantage  should 
be  reversed  for  two  and  three-dimensional  problems. 

The  second  problem  considered  is  the  forced  vibration  of  a 
nonlinear  string.  In  this  case,  the  exact  solution  is  chosen  to  be  a 
product  of  a  spatial  term  and  a  temporal  factor: 


Z|  *  [i/(n+1  )]*sin(t) 
x,  *  [(n-i)/(n+1 )]*sin(t) 


for  i  <s  (n+1)/2 
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for  i  >  (n+1 )/2 

This  displacement  pattern  is  a  symmetric  bilinear  function  of  the 
spatial  coordinate,  vith  zero  displacement  at  the  ends  of  the  string, 
and  a  marl mum  value  of  1/2  at  the  center  of  the  string.  Again,  the 
string  was  subdivided  into  fifty  point  masses,  and  six  modes  were  used 
in  the  reduced  algorithms. 

The  nonlinearity  present  in  this  problem  is  a  simple  quadratic 
nonlinearity  on  the  diagonal  of  the  stiffness  matrix: 

K  a  tridiag(-1 ,  2«(1+x?),  -1) 


Qiven  this  construction  of  the  exact  solution,  the  accuracy  of  the 
numerical  schemes  can  be  evaluated.  The  analysis  was  performed  with 
three  different  time  steps,  and  the  errors  are  presented  in  Table  2. 
Displacement  errors  associated  with  the  shortest  time  steps  are 
displayed  in  Figure  10.  Note  that  the  most  serious  errors  are  present 
at  the  center  of  the  string:  this  is  due  to  the  "kink"  in  the  exact 
solution  at  the  this  point.  This  disturbance  propagates  along  the 
charaoteristics  in  the  fully  continuous  problem,  but  is  "washed-out" 
by  a  numerical  procedure  (such  as  those  used  in  this  example)  that 
ignores  the  presence  of  characteristic  ourves.  (This  phenomenon  is 
similar  to  what  occurs  when  attempting  to  solve  hyperbolio  partial 
differential  equations  using  finite-difference  methods. )  Away  from  the 
oenter  of  the  string,  the  solutions  are  quite  accurate. 
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A  simple  nonlinear  oorreotion  scheme  of  subtracting  the  residual  at 
the  l**1  step  from  the  load  vector  applied  at  the  (i+1)st  step  was  used 
with  considerable  success  in  this  problem.  Further  iteration  within 
the  time  step  decreased  the  error  for  the  direct  method  (refer  to 
Table  2),  but  did  not  Improve  the  reduced  algorithms. 


The  last  example  presented  is  the  problem  of  free  vibration  of  a 
nonlinear  string.  The  initial  conditions  are  the  same  as  in  the  last 
problem,  but  the  nonlinearity  is  weaker: 


K  >  tridiag(-1,  2+xj*,  .1) 


The  results  of  the  three  integration  schemes  are  shown  in  Figures 
11,  12  and  13.  Maximum  values  of  the  displacement,  velocity,  and 
acceleration  are  given  in  Table  3,  along  with  more  refined  values 
obtained  by  dlreat  integration  using  a  small  time  step  (&t  =  0.01). 


Til.  Conclusions  and  Recommendations  for  Further  Study 


This  researoh  is  intended  only  as  a  first  step  towards  a  general 
algorithm  for  the  solution  of  time-dependent  nonlinear  problems  of  the 
form  given  by  Eqn  9.  Considerable  work  needs  to  be  done  in  a  variety 
of  areas  in  order  to  achieve  a  robust,  efficient  alternative  to  the 
direct  integration  of  the  full  equations  of  motion.  The  example 
problems  presented  are  limited  in  scope  and  geometry,  and  a  number  of 
two  and  three-dimensional  problems  involving  nonlinear  geometrio 
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and/or  material  behavior  are  being  considered.  These  more  realistic 
problems,  and  the  topics  considered  next,  form  the  outline  of  the 
Ph.D.  research  of  the  first  author. 

There  are  four  areas  in  which  further  research  and  development  is 
presently  being  concentrated: 

1 )  Development  of  efficient  and  reliable  iteration  schemes 

The  numerical  integration  of  the  nonlinear  equations  of  motion 
across  a  time  step  is  equivalent  to  the  solution  of  a  system  of 
nonlinear  equilibrium  equations  at  the  end  of  the  time  step.  This 
solution  is  obtained  using  some  sort  of  iteration  strategy,  involving 
formation  of  the  equations  of  motion  at  the  end  of  the  step,  and 
checking  for  a  residual  error.  Because  one  of  the  most  expensive 
operations  involved  in  the  proposed  algorithm  is  the  formation  of  the 
reduced  stiffness  at  each  Iteration,  the  number  of  iterations  required 
to  satisfy  the  equilibrium  equations  must  be  kept  to  a  minimum  in 
order  to  achieve  the  highest  computational  efficiency.  This  requires 
that  the  iteration  scheme  have  rapid  convergence  characteristics. 
Present  research  is  being  focused  on  the  Newton  family  of  methods  as 
an  alternative  to  a  iterative  application  of  Newmark's  Method. 
Particular  attention  will  be  given  to  the  use  of  the  algorithm  with 
comprehensive  material  inelastioity  models,  such  as  the  bounding 
surfaoe  model  for  soils  [11].  Means  will  be  sought  for  combining  the 
inherent  need  of  suoh  models  for  iteration  and  integration  across  a 
given  time  step  with  the  Iteration  requirement  of  the  dynamic 
analysis. 


2)  Increased  efficiency  in  formation  of  reduced  problem 

The  limiting  factor  in  the  present  implementation  of  the  reduced 
algorithm  is  the  process  of  forming  the  reduced  matrices.  There  are  a 
number  of  avenues  open  for  refinement  of  this  process,  including  an 
asymptotic  approach  (see  [10]),  the  construction  of  the  reduced 
matrices  at  the  element  level,  and  the  use  of  parallel  processing. 
Progress  on  this  topic,  coupled  with  an  efficient  iteration  scheme  as 
mentioned  in  the  last  section,  would  greatly  speed  the  proposed 
algorithm. 

3)  Adaptive  determination  of  number  of  modes  required 

The  appropriate  number  of  modes  required  for  an  acceptable  solution 
may  not  be  apparent  to  the  analyst  either  before  or  during  the 
solution  process.  In  addition,  the  determination  of  whether  the 
reduced  coordinates  need  to  be  updated  must  be  automated,  in  order  to 
be  able  to  use  the  proposed  algorithm  in  any  way  resembling  the  use  of 
a  ’blaok  box'.  Both  these  issues  need  to  be  addressed  before  the 
reduced  algorithm  can  compete  with  direct  integration  schemes,  and 
both  require  further  research  into  the  effect  of  projection  error  in 
the  reduction  prooess. 

4)  Modifications  for  stiff  (i.e.  soil-structure)  systems 

Preliminary  results  from  this  research  indicate  considerable 
promise  for  the  reduced  algorithm  for  solving  certain  stiff  problems. 
In  particular,  soil-structure  interaction  problems  are  characterized 
by  such  a  wide  range  of  frequencies  that  a  representative  time  scale 


for  the  analysis  of  the  coupling  between  the  soil  and  structural 
components  is  difficult  or  impossible  to  find.  One  promising 
possibility  would  involve  projecting  the  soil-structure  problem  onto 
the  soil  deformation  modes  that  impart  energy  to  the  structure,  and 
the  lowest  modes  of  vibration  of  the  structure  itself.  The  Lanczos 
algorithm  has  been  used  by  the  authors  with  an  origin  shift  to  obtain 
reduced  coordinates  for  such  problems,  and  work  is  underway  on 
incorporating  these  coordinates  into  the  proposed  algorithm. 


In  conclusion,  the  use  of  a  reduced  set  of  modal  coordinates  offers 
an  alternative  to  direct  integration  of  the  full  equations  of  motion. 
The  research  presented  in  this  work  represents  a  first  step  towards  an 
efficient  implementation  of  this  reduction  for  a  variety  of  nonlinear 


problems 
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Appendix  1 :  Table  of  Symbols 


Symbol 


Definition 


b 

D 

f 

g 

b 

K 

M 

n 

P 

Qp 

R 

t 

Up 

wi 

x 

X 

X 

z 


Average  bandwidth  of  stiffness  matrix 
Diagonal  matrix  of  squared  frequencies 
Finite  element  load  vector 
Reduced  load  vector 
Component  of  f  orthogonal  to  g 
Finite  element  stiffness  matrix 
Finite  element  mass  matrix 
Number  of  equilibrium  equations 
Number  of  reduced  coordinates 
Individual  trial  mode  shape 
Orthogonal  matrix  of  trial  mode  shapes 
Reduced  stiffness  matrix 
time  variable 

Orthogonal  matrix  of  exact  eigenvectors 

Modal  natural  frequency 

Vector  of  nodal  displacements 

Vector  of  nodal  velocities 

Vector  of  nodal  accelerations 

Vector  of  reduced  coordinates 


z 


Veotor  of  reduced  accelerations 


Table  1 :  Errors  in  Solution  of  Linear  Free  Vibration  Problem 


Type  of 

Analysis 

Maximum  Absolute 
Displacement  Velocity 

Errors 

Acceleration 

Cost 

(cpu-sec) 

Direct  integration 

0.0424 

0.0203 

0.0139 

9.53 

6  Eigenvectors 

0.0433 

0.0173 

0.0126 

35.70 

10  Eigenvectors 

0.0426 

0.0210 

0.0145 

80.70 

6  Lanczos  vectors 

0.0413 

0.0202 

0.0137 

32.63 

10  Lanczos  vectors 

0.0465 

0.0169 

0.0118 

77.14 

Table  2:  Errors  in  Solution  of 

Nonlinear 

Forced  Vibration  Problem 

Type  of 

Maximum 

Absolute 

Errors 

Cost 

Analysis  Displacement 

Velocity 

Acceleration 

(cpu-sec) 

Direct,  128  steps 

0.0655 

0.0860 

0.1092 

18.71 

Direct,  64  steps 

0.1014 

0.1411 

0.1894 

9.85 

Direct,  32  steps 

0.2101 

0.1749 

0.1837 

5.57 

Direct  w/  iteration, 128  steps 

0.0098 

0.0107 

0.0128 

25.74 

Direct  w/  iteration,  64  steps 

0.0253 

0.0277 

0.0321 

13.59 

Direct  w/  iteration,  32  steps 

0.0942 

0.0964 

0.1083 

7.90 

6  Eigenvec.,  128  steps 

0.0440 

0.0434 

0.0698 

65.61 

6  Eigenvec. ,  64  steps 

0.0751 

0.0673 

0.1243 

34.71 

6  Eigenvec. ,  32  steps 

0.2016 

0.1975 

0.1936 

19.62 

6  Lanczos  vec.,  128  steps 

0.0371 

0.0313 

0.0542 

63.42 

6  Lanczos  vec.,  128  steps 

0.0644 

0.0657 

0.0961 

31.97 

6  Lanczos  vec.,  128  steps 

0.2306 

0.2015 

0.2049 

16.82 

Table  3:  Approximate  Solution  of  Nonlinear  Free  Vibration  Problem 

Cost 

(cpu-sec) 

Type  of 

Analysis 

Maximum  Absolute 
Displacement  Velocity 

Response 

Acceleration 

Direct  integration 

0.8440 

0.5135 

0.5982 

12.94 

Direct  v/  small  step 

0.8177 

0.4902 

0.5521 

81.21 

6  Eigenvectors 

0.8572 

0.5601 

0.6674 

60.41 

6  Lanczos  vectors 

0.8650 

0.5651 

0.6762 

57.79 

(All  computation  times  are  for  a  VAX  11/750  with  floating-point  hardware) 
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Exact  eigenvectors 
for  example  1 


(a)  Modes  1  -  4 


(b)  Modes  5-7 


(o)  Modes  8-10 


Figure  2 


Lanczos  vectors 
for  example  1 

(a)  Modes  1  -  4 


(b)  Modes  5-7 


(o)  Modes  8-10 


Figure  3 

One-dimensional  solution  history 


Displacement 

Velocity 

Acceleration 


(b)  Geometry  of 

solution  display 


Numerical  solution  history 
for  the  linear 
free-vibration  problem 
(Direct  integration) 

(a)  Approximate  displacement 


Figure  6 


Numerical  solution  history 
for  the  linear 
free-vibration  problem 
(6  Eigenvectors) 


(a)  Approximate  displacement 


MxiNM  absolute  disf lacomnt  :  6.372 


HwiMUM  absolute  velocity  :  0.4902 


(b)  Approximate  velocity 


(o)  Approximate  acceleration 


Maxi  mum  absolute  acceleration  :  0.02606 


Figure  7 
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Numerical  solution  history 
for  the  linear 
free-vibration  problem 
(6  Lanczos  vectors) 


(a)  Approximate  displacement 


► 


(b)  Approximate  velocity 


f 


(o)  Approximate  acceleration 


Figure  8 


Displacement  errors 
for  the  linear 
free-vibration  problem 


(a)  Direct  integration 


(b)  6  Eigenvectors 


(o)  6  Lanozos  Vectors 


Figure  9 


Figure  10 


Displacement  errors 
for  the  nonlinear 
forced  vibration  problem 
(example  2) 

(a)  Direct  integration 


(b)  6  Eigenvectors 

i 


(o)  6  Lanczos  Vectors 


Figure  11 


Numerical  solution  history 
for  the  nonlinear 
free-vibration  problem 
(Direct  integration) 


(a)  Approximate  displacement 


(b)  Approximate  velocity 


(c)  Approximate  acceleration 


Figure  12 


Figure  13 


Numerical  solution  history 
for  the  nonlinear 
free-vibration  problem 
(6  Lanczos  vectors) 


(a)  Approximate  displacement 


